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RELAXATION IN A-BODY SIMULATIONS OF DISK GALAXIES 
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ABSTRACT 

I use A-body simulations with two mass species of particles to demonstrate that disk galaxy simulations 
are subject to coUisional relaxation at a higher rate than is widely assumed. Relaxation affects the vertical 
thickness of the disk most strongly, and drives the velocity ellipsoid to a moderately flattened shape 
similar to that observed for disk stars in the solar neighborhood. The velocity ellipsoid in simulations 
with small numbers of particles quickly approaches this shape, but shot noise also dominates the in- 
plane behavior. Simulations with higher, but reachable, numbers of particles relax slowly enough to be 
considered collisionless, allowing the in-plane dispersions to rise due to spiral activity without heating the 
vertical motions. Relaxation may have affected many previously published simulations of the formation 
and evolution of galaxy disks. 

Subject headings: Galaxies: kinematics and dynamics — numerical methods 



1. INTRODUCTION 

Simulations over several decades have contributed a great 
deal to our understanding of the complex dynamical be- 
havior of stellar systems. The numbers of particles em- 
ployed in typical work has risen steadily with the avail- 
able computational power, but even today we are unable 
to employ as many particles as there are stars in a galaxy, 
at least on a routine basis. Thus particles in simulations 
are typically fewer in number and more massive than are 
stars in a galaxy, and therefore the density distribution 
suffers from a higher level of shot noise than expected in 
t he stellar component o f the system being mimicked. 

IChandrasekharl (|1941[ ) estimated the relaxation time, or 
the time required for deflections by other stars to cause a 
star to lose all memory of its original trajectory, to be 

where fi is the mass of a star that moves at the typical 
speed V and n is the number of stars per unit volume. The 
ratio A ~ Re/bmim with Re being the half-mass radius of 
the galaxy, and ^min 2G/i/u^, the impact paramater 
below which the approximation that deflections are small 
fails badly. The Coulomb logarithm. In A, indicates that 
every decade of increase of the impact parameter makes 
an equal contribution to the relaxation rate, because the 
number of stars rises with distance to compensate for the 
diminishing influence of each individual encounter. In sim- 
ulations, one generally replaces b^i^ with the gravity soft- 
ening length, e, but distant encounters are unaffected by 
softening and fully contribute to relaxation. The reason 
for softening the interparticle force is mostly to avoid the 
need for shorte r time steps during close approaches be- 
tween particles (jBinnev fc Trem ainc 2008. hereafter BT08, 
pl23). 

Following BT08, we set ~ GNfi/Re, n ~ iN/iAnRl), 
and define a dynamical time Tdyn = Re/ v. With these 
substitutions, eq. ([T]) becomes 
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indicating that star-star encounters are utterly negligible 
for galaxies. The effect of a halo of particles whose masses 



are much smaller than those of stars is further to lengthen 
the relaxation time, since it increases v without contibut- 
ing to the deflections. Eq. ^ also suggests that simula- 
tions with N > 10^ should be safely collisionless, especially 
as the In N term in the denominator is an overestimate for 
the values of e usually employed. 

While scattering in disks is more rapid, see §2, relaxation 
caused by star-star encounters remains much too slow to 
determine the random speeds of stars in the solar neigh- 
borhood. The local dis t ribution of peculiar ve locities {e.g. 
iNordstrdm et al\ |2004 iHolmberg eTaU l2009t) has a flat- 
tened tri-axial shape. The dispersion in the direction of 
the Galactic center, a^, is the largest component, with 
that in the orbital direction, ct^, being about 70% as great 
as expected from epicycle theory (BT08, pl70), while the 
dispersion, a^, normal to the Galactic plane is some 60% 
of au , and is the smallest; 

In an important paper, llda et al\ (11993') showed that the 
shape of the velocity ellipsoid is consistent with scatter- 
ing by dense mass clumps, such as giant molecular clouds 
(GMCs). Their predicted flattening depends on the slope 
of the Galactic rotation curve, but were this locally flat, 
they predict cr^ ~ O.Gctr, as observed. While GMCs read- 
ily redirect peculiar motions, they are less efficient at in- 
creasing them; G MCs are believed to be insufficiently m as- 
sive or numerous (lLacevlll991l: lH annmen fc Flvnn|[2002l ) to 
create the high dispersion of the oldest thin-disk stars. 

Disk stars are also believed to be scattered at the Lind- 
blad resonances of transient spirals (jCarlberg fc SellwoodI 
|1985() . which increase only the in-plane dispersions. The 
vertical dispersion should be unaffected because the rapid 
vertical oscillation of a star is adiabatically invariant at the 
slow r ate at which it encounters the spiral density varia- 
tions (jCarlberd 119871) . Thus spiral scattering is needed 
to account for the high peculiar speeds of the older thin- 
disk stars while massive gas clumps redirect the peculiar 
velocities to ma intain the shape of the velocity ellipsoid 
(|Sellwoodll20Tl . 

Most simulations of galaxy disks do not include a sepa- 
rate population of heavy particles, and therefore have no 
agent to redirect the random motions generated by spiral 
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waves into the vertical direction if they we re truly col- 
lisionl ess. However, a few authors, n otably iQuinn et al\ 
(|1993[ ) and lMcMillan fc DehnenI (|200iD . have worried that 
isolated disks thicken in simulations as spiral activity heats 
the in-plane motions. Here I show that coUisional relax- 
ation is the likely origin of their finding, and should be a 
concern for all simulations of disk galaxies. 

2. RELAXATION IN DISKS 

Formulae ([T]) and ([2]) were derived assuming a mass dis- 
tribution that is roughly uniform in 3D. The relaxation 
rate in disks dif fers for a t least four reasons, as were mostly 
pointed out bv lRvbickil (|1972D . 

First, because disks are rotationally supported, stars 
pass each other at speeds that are a small fraction of the 
circular orbit speed, say with /3 ~ 0.1. The predomi- 
nantly slow encouters mean that the typical deflection rate 
is boosted by the factor and the time needed for the 
random impulses to amount to a peculiar speed of /3v is 
shorter by a factor /3^. 

Second, eq. ^ results from an integration of mean- 
square deflections over all impact parameters. The 3D 
integration volume element becomes an element of area 
in 2D, and therefore the impact parameter enters with 
one power less. Thus the Coulomb logarithm is replaced 
by the factor — &max); indicating that scattering is 

dominated by close encounters. Distant encounters are 
negligible in disks because the number of stars does not 
rise rapidly enough with distance to compensate for the 
diminished force each exerts. 

Real galaxy disks are neither razor thin, nor spherical. 
In this case, the spherical dependence applies at ranges up 
to the typical disk thickness, zq, beyond which the contri- 
bution to scattering drops quickly. Thus we should replace 
Re in the Coulomb logarithm by zq, which slightly reduces 
the relaxation rate. 

Third, the local number density of stars is higher so 
that N ~ ttRIzqu, which increases the density to be used 
in eq. ([Ij by the factor Rs/zq. This factor causes another 
order of magnitude increase in the relaxation rate. 

Combining these three considerations yields a relaxation 
time in disks that is shorter than in 3D distributions by 
the factor 

3 / Zq \ In {Re/bnun) , . 

\Re) ln(2o/Vin)' ^ ' 

which is typically a few xlO^**! It should be noted that 
the cube of /? arises from the time needed for scattering 
to produce a peculiar velocity i.e. equal to that of a 
typical disk star; the relaxation rate is increased over the 
3D rate by only a single factor of (3, but together with the 
higher density, the increase is still ^ 100-fold. 

The fourth factor in real disks is the existence of GMCs, 
whose role in determining the shape of velocity ellipsoid 
was described in the introduction. 

3. SIMULATIONS 

Here I report a few simulations to test these theoretical 
predictions; a fuller study will be presented elsewhere. 

Since simulations of all realistic disks heat due to insta- 
bilities, simple measurements of the heating rate will not 
det ermine the relaxation rate directly. Therefore, follow- 
ing |Hoh3 (|l973i ). I employ separate particle species having 



different masses, since energy exchange between particles 
of differing mass is a reliable indicator of relaxation. 

The smooth, half-mass Mestel disk (BT08, p99) was 
proved bv.Toomre ( 1981) to have no global instabilities, al- 
though [SellwoodI (f2012[ ) found subtle instabilities, caused 
by non-linear effects, that gave rise to spiral activity in 
large- TV models. It seems unlikely that 3D motion would 
alter the unusual stability properties of this disk, which 
make it more attractive for this test than other more generic 
disks that could host global linear instabilities. 

3.1. Setting up the disk 

In order to construct an equilibrium model, I start from 
the two-integral DF for the razor-thin disk (jZand [19761 : 

lToomrd[T977[) : fz{E,L^) cx Lle'^/^R, where £; is a par- 
ticle's specific energy and its specific z-angular mo- 
mentum. The parameter q determines the radial velocity 
dispersion = Vo{l + q)~^^'^ , with Vb being the circu- 
lar orbital speed at all radii. The value of Toomre's local 
stability parameter for a razor-thin Mestel disk is 



^~^R.min 3.36/(1 + g)l/2' 

where / is the active mass fraction in the disk; both and 
Q are independent of radius. In order obtain a disk with 
/3 < 0.15, which must also have Q > 1, I adopt / = 0.25, 
with the remaining mass in a rigid halo. 

I apply inner and outer tapers to limit the radial extent 
of the disk: foiE,L,) = fz[l+iL^/ +{L J L,)^]-\ 
where Li and Lg are the central angular momentum val- 
ues of the inner and outer tapers respectively. I choose 
Lo — 15Li, and further restrict the extent of the disk by 
eliminating all particles whose orbits would take them be- 
yond 20i?i, where Ri = Li/Vo is the central radius of the 
inner taper. With these tapers Re — 8Ri. 

I thicken the disk by giving it the Gaussian vertical 
density profile p{R,z) = cxp(— z^/2zQ)/(27rzo), with 
being the vertically integrated surface mass density 
at radius R. I estimate the equilibrium vertical velocity 
dispersion at each z-height by integrating the ID Jeans 
equation (BT08, eq. 4.271) in the numerically-determined 
potential. I adopt an unrealistically small value zq = 
0.05i?i, that is independent of R, in order to reveal the 
effects of relaxation clearly. 

The active mass gives rise to a weaker central attraction 
in the midplane than that which should arise from a razor- 
thin, infinite, full-mass Mestel disk. In order to maintain 
equilibrium, I add a rigid central attraction to the grid- 
determined forces from the particles in the disk at each 
step. This unchanging, spherically symmetric, rigid cen- 
tral attraction is pre-tabulated from differencing —Vq/R 
from the grid-determined attraction of a smooth density 
created from the thickened, tapered disk. 

3.2. Numerical details 

I determine the gravitationa l field of the particles usin g 
the 3D polar grid described in lSellwood fc Valluril ()1997[ ). 
The grid has 200 rings, 256 spokes, and 375 vertical planes, 
and I adopt a cubic spline softening rule that yields the full 
attraction of a point mass at distances > 2e. I choose Ri = 
5 grid units, the grid planes are 0.027?^ apart, e = 0.025i?i, 
a basic timestep of 0.025i?i/Vb, and advance particles at 
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Fig. 1. — Measurements over a broad swath of the disk around 
the half-mass radius of the time evolution of the average an (left) 
and CTz (right). The top row shows results with N = 2x 10* in each 
species and the particle number was increased 10-fold from row to 
row. The red (green) curves apply to the heavy (light) particles 
respectively. Velocities are in units of Vq, distances in units of Ri, 
and the rotation period at Re is ~ 50 in these time units, where 
G = Vo = Ri = l. 



radii R > 2Ri at intervals that double in duration from 
this radius and again with every factor 2 increase in R. 
The duration of the simulations was ~ 20 full rotations of 
the disk at i?e , or ^ 3 Gyr when scaled to the Milky Way. 

In order to create two populations of particles with un- 
equal masses, I employ every particle selected from the DF 
twice, placing the two particles each at separately cho- 
sen random azimuths, and make one 9 times more mas- 
sive than the other. The total masses of both populations 
are set to yield the desired total disk mass (/ = 0.25 in 
eq. SI) and a combined Q — 1.5, which corresponds to 
(jR ~ 0.14Vo and cr^ ~ O.IVq, i.e. in the ratio expected for 
a flat rotation curve (BT08). Thus fitotai/Vb = ^ ~ 0.17. 

3.3. Results 

Figure [1] shows the evolution of (Th and , averaged over 
a broad radial range centered on Re, measured for both 
the heavy (red) and light (green) particles in three simula- 
tions. Each row is from a separate simulation, differing in 
the number of particles, as indicated in the left-hand pan- 
els. Since the simulations start with randomly placed par- 
ticles, the i nitial behavior is dominat ed by swing-amplified 
shot noise ijToomre fc Kalnaiill99lh . and I stop the cal- 
culations after ~ 20 rotations at Re, while spiral activity 
is on-going. 

The radial velocity dispersion rises due to spiral activity, 
which has lower initial amplitude as the number of par- 
ticles rises. The vertical dispersion and rms z-thickness 
(Figure [2]) of each population rises rapidly in the smallest 
simulation, but are almost constant in the largest, demon- 



FlG. 2. — Continuation of Fig. [T] rms ^-thickness (left) and the 
ratio Uz/crji (right). 



strating that the models were in initial equilibrium. How- 
ever, the tendency for mass segregation, unmistakable in 
the top row, is still visible in the bottom row. 

The different heating rates for the two mass species is 
clear evidence for relaxation, as is also the rapid rise in the 
vertical motions of both populations as the number of par- 
ticles is decreased. Notice also, from the top right panel of 
Figure [5] that the velocity ellipsoid shapes of both popula- 
tions rapidly become rounder. The larger experiments, on 
the other hand, manifest a slower energy equipartition rate 
and also support the theoretical prediction that spirals do 
not heat vertical motions. 

The rate at which rises due to encounters, one step 
before arriving at eq. ([T}, is 

dt V 

In the simulations, ~ Md/Nh with being the num- 
ber of heavy particles, A = zg/e, v = /3Vb, and n 
Nh/{TrRlzo). We set Ma = fMtot, and Vq^ = GMtot/{Re), 
and assume that the vertical velocity dispersion, az , is ris- 
ing due to relaxation. Thus we expect 

dt 



(5) 



40 



(6) 



NhPzo 

if = 1, / = 0.25, zo = i?«/20 = 2e, and P ~ 0.17. 

I estimate the initial value of da'^/dt 1.1 x 10^^ from 
the light particles in the intermediate simulation. This is 
to be compared with 40/7V,i = 2 x 10"'', since 7V/i = 2x10^. 
Thus the predicted scattering rate is about 20 times that 
observed. This discrepancy appears to be due to some 
additional grid smoothing in the simulation. Relaxation 
is even slower when softening and grid cell sizes are in- 
creased, and also when fewer sectoral harmonics, to, con- 
tribute to the grid-determined forces. In Figs. [T] & [2J the 
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force calculation discarded contributions to the forces from 
m > JTiinax = 32, and the relaxation rate in further simu- 
lations drops by about 30% for each factor 2 reduction in 
"T-max- Applying a cut-off in the expansion of the density 
distribution effectively smooths each particle in the az- 
imuthal direction, weakening the attraction between disk 
particles and slowing the relaxation rate. Since relaxation 
is inhibited by the extra smoothing in grid methods, I 
have deliberately employed a finer grid and more sectoral 
harmonics than is normally required in order to highlight 
relaxation. 

A dependence on grid smoothing is physically reason- 
able, b ut contrasts with the finding of lHernguist fc BarnesI 
(|1990( ) that the relaxation rate does not depend on the 
simulation method. Note that they were testing spherical 
models, where relaxation is dominated by density fluctua- 
tions on large scales and is little unaffected by short-range 
smoothing, which is the reason they found the rate of re- 
laxation to be the same in every valid numerical method. 
However, relaxation in disks is dominated by close encoun- 
ters (§2), making the rate much more dependent on the 
degree of smoothing. Therefore, it seems reasonable to 
expect a higher relaxation rate for the same N in direct 
methods, e.g. tree codes, than in grid-based methods. 

4. DISCUSSION AND CONCLUSIONS 

Both the theoretical argument in §2 and the numerical 
results in §3 confirm that two-body scattering in disks is 
much more rapid than that expected in a spheroidal model 
with the same N. This is because disk particles pass each 
other at speeds much lower than the orbital speed, and 
because the density of particles in a disk is higher than 
were the same number spherically distributed. Further- 
more, unlike in 3D systems, relaxation is dominated by 
the particles within a few disk thicknesses, and the cumu- 
lative effect of distant encounters is less important. 

The experiments also confirm that spirals, which heat 
the in-plane motions (left panel of Fig. [T]), do not cause 
even extremely thin disks to thicken, when relaxation is 
slow enough (bottom row). With fewer particles, relax- 
ation causes the shape of the velocity ellipsoid of both 
species to evolve rapidly to roughly the shape observed in 
t he solar neighborhood (top r ight panel of Fig. [5]). 

[McMillan fc DehnenI ()2007t ) tried to investigate why the 
disks in their simulations thickened. They found that 
thickening was inhibited in parallel experiments in which 
the position and velocities of each disk particle were ro- 
tated through a random angle after each time step. Having 
suppressed collective spiral responses by this stratagem, 
they concluded that thickening was due to spiral heat- 
ing. However, two-body relaxation arises from the time- 
integrated perturbing forces as particles pass, which must 
also have been suppressed in this test. 

Possible relaxation needs to be taken into account when 
interpreting results from simulations in a wide variety of 
contexts. For example, the number of star particles gen- 
rally employed in disk formation simulations is in the range 
of a few ^ 10^, a s reported in the code comparison by 
IScannapieco et all ((2012) . Simulations of star-forming disks 
could be particularly severely affected by relaxation, since 
new stars fo r med f rom cold gas will have the smallest /?. 
iHouse et al\ (|201ll ) compare the thickening of their sim- 



ulated disk with SDSS data, but it is unclear what this 
test shows because relaxation probably caused some thick- 
ening in their model. Relaxation may have contributed 
to disk thick ening in the simu l ation of dwarf disk galaxy 
formation bv iGovernato et all (|2010D . which employed 
5 X 10^ particles. Other areas where relaxation may have 
affected the co nclusions include studies of the survival o f 
thin disks (e.g. 'Robert son et aLl [20061 : iMoster et aLll2010l ) 
and thickness variations during radial migrat i on in disks 
(e.g. iMinchev et al\ 120111: iRoskar et all 120121: iBird et al\ 
l2013f ). In partic ular, the actions of particles, when calcu- 
lated exactly bv lSolwav et al\ (|2012D . were conserved only 
on average in their simulations with N > 10^; it seems 
likely that the inexact conservation of this quantity could 
have been due to slow relaxation. 

Since disk thickening is most strongly affected by re- 
laxation, the outcome of a simulation will depend on the 
number of particles employed. On the one hand, the shape 
of the velocity ellipsoid is determined by relaxation in sim- 
ulations with modest N. The heavy particles scatter each 
other somewhat as GMCs scatter stars in galaxies so that, 
paradoxically, low-quality simulations get the right shape 
of the velocity ellipsoid for the wrong reason! But in small- 
N simulations, the in-plane dynamics is dominated by the 
collective responses to shot noise. 

With large N, on the other hand, coUisionless in-plane 
dynamics is more faithfully represented and coupling of 
spiral heating to the vertical motion is weak. But to mimic 
the evolution of the velocity ellipsoid, one would need to in- 
clude a population of extra-heavy particles to redirect the 
in-planc motions. Note that the se heavy particles would 
also affect the in-plan e dynamics (jToomre fc Kalnaisll99ll : 
iD'Onghia eta!\mi^ . 

Thus any simulation of an isolated stellar disk that thick- 
ens probably does so through 2-body relaxation. If gas is 
included, clumps of gas particles may behave as scattering 
centers, as do the GMCs in real disks, but with a smaller 
mass ratio to the star particles. Simulations that mimic 
the full hierarchical evolution will have other sources of 
vertical heating, such as in-falling dwarf galaxies and sub- 
halos, etc. However, segregation of star particles of differ- 
ent masses would remain a valid diagnostic of relaxation 
in all these contexts. 
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